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We  study  the  two  dimensional  Kolmogorov  flows  in  the  limit  as  the  forcing  frequency 
goes  to  infinity.  Direct  numerical  simulation  indicates  that  in  this  limit  the  low  frequency 
energy  spectrum  evolves  to  a  universal  A:-4  decay  law.  We  derive  effective  equations  gov¬ 
erning  the  behavior  of  the  large  scale  flow  quantities.  We  then  present  numerical  evidence 
that  with  smooth  initial  data,  the  solution  to  the  effective  equation  develops  a  k~4  type 
singularity  at  a  finite  time.  This  gives  a  convenient  explanation  for  the  k~4  decay  law 

exhibited  by  the  original  Kolmogorov  flows. 
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§  1.  Introduction 


Toward  the  end  of  the  50’s,  Kolmogorov  proposed  to  study  the  two-dimensional  incompress¬ 
ible  flow  with  a  high  frequency  forcing  term  as  a  model  to  help  understand  the  transition  to 
turbulence  and  the  inverse  cascade  process  in  a  turbulent  flow.  These  so-called  Kolmogorov 
flows  are  described  by  the  following  equation: 


(1.1) 


ut  4-  u  •  Vu  +  Vp  =  vAu  +  uk2j 
V  •  u  =  0. 


sin(fc/a:2) 

0 


Here  v  is  the  dynamic  viscosity,  kj  is  the  forcing  frequency.  Since  then,  this  problem  has 
attracted  a  lot  of  attention.  Considerable  progress  has  been  made  on  both  the  experimental 
and  theoretical  side.  We  refer  to  [10]  for  a  discussion  on  the  laboratory  realization  of 
the  Kolmogorov  flows.  The  theoretical  side  of  this  problem  has  mostly  been  concentrated 
on  the  hydrodynamic  stability  viewpoint.  Linear  stability  of  the  basic  Kolmogorov  flow 
(a  stationary  ear  flow  proportional  to  the  forcing)  was  studied  in  [6].  The  amplitude 
equations  for  the  weakly  nonlinear  theory  were  derived  in  [9,  15].  It  was  found  that  this 
problem  belongs  to  a  wider  class  of  problems  which  exhibit  large  scale  instability.  Later  stage 
development  of  the  Kolmogorov  flow  was  extensively  studied  more  recently  in  [11,  12]  using 
a  dynamical  system  viewpoint.  Several  scenarios  of  the  transition  to  chaos  and  turbulence 
were  identified  by  resorting  to  careful  numerical  simulations. 

In  this  paper,  we  take  a  different  approach  and  study  Kolmogorov’s  problem  in  the 
following  asymptotic  limit: 


Our  main  interest  is  to  understand  the  inverse  cascade  process  in  the  Kolmogorov  flow.  Our 
extensive  numerical  results  indicate  that,  even  if  we  start  with  a  smooth  initial  data  with 
an  exponentially  decaying  energy  spectrum,  the  high  frequency  forcing  term  pumps  energy 
into  lower  frequencies  so  that,  beginning  at  a  finite  time,  the  energy  spectrum  exhibits  a 
universal  k~A  algebraic  decay  law  at  the  low  frequencies. 
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To  explain  this  phenomenon,  we  derive  the  effective  equations  governing  the  evolution  of 
the  large  scale  flow  quantities.  We  present  strong  numerical  evidence  that  the  solution  of  the 
effective  equation  with  smooth  initial  data  form  a  singularity  at  a  finite  time.  Furthermore, 
the  energy  spectrum  at  the  time  of  singularity  formation  is  precisely  k~4. 

We  remark  that  in  this  paper  the  term  “cascade”  is  used  loosely  to  describe  situations 
where  either  energy  or  enstrophy  is  shifted  in  the  frequency  space.  We  do  not  have  evidence 
that  energy  or  enstrophy  is  transferred  continuously  in  the  frequency  space  with  nonzero 
flux  at  each  wavenumber  since  we  did  not  measure  these  quantities. 

This  -\per  is  organized  as  follows.  In  the  next  section  we  present  results  of  direct 
numerical  simulation  of  the  Kolmogorov  flow  in  the  limit  as  the  forcing  frequency  goes 
to  infinity.  The  emphasis  will  be  placed  on  the  evolution  of  the  energy  spectrum.  In 
Section  3,  we  derive  the  effecetive  equations  for  the  averaged  quantities  using  multiple  scale 
asymptotics.  In  Section  4,  we  present  our  numerical  evidence  for  the  singularity  formation  in 
the  solutions  of  the  effective  equations.  We  conclude  this  paper  in  Section  5  by  commenting 
on  the  relevance  of  our  findings  for  more  general  flows,  including  three-dimensional  flows. 

§  2.  The  Inverse  Cascade  Phenomenon  in  the  Kolmogorov  Flows 


As  we  mentioned  in  the  introduction,  we  are  interested  in  studying  the  inverse  cascade 
phenomena  associated  with  the  Kolmogorov  flow  in  the  following  asymptotic  limit: 


(2.1) 


v  = 


1 

kf 


=  e. 


As  can  be  easily  seen  from  the  derivations  presented  in  the  next  section  (see  Remark  3.2), 
this  is  the  only  distinguished  limit  for  which  both  the  viscous  force  and  the  external  high 
frequency  force  play  important  roles  in  the  effective  equations  governing  the  evolution  of 
large  scale  structures. 

Under  the  scaling  (2.1),  equation  (1.1)  takes  the  following  form: 


(2.2) 


u\  4-  ue  •  Vu1  +  =  sAu£  +  i 

V  ■  ue  =  0 


(  sinf  \ 

V  0  / 
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To  study  the  influence  of  the  small-scale  forcing  on  the  evolution  of  the  large  scales,  in 
particular,  to  study  the  evolution  of  the  energy  spectrum  associated  with  the  Kolmogorov 
flow,  we  have  numerically  integrated  (2.2)  using  a  number  of  initial  velocity  fields  and 
studied  the  limit  as  £  — ►  0.  We  discretize  (2.2)  using  a  Fourier-collocation  method  in  space, 
and  a  third  order  Runge-Kutta  method  in  time.  More  details  of  the  numerical  method  will 
be  given  in  Section  4.  Our  calculation  differs  from  the  previous  calculations  of  Sulem  et.al. 
[16]  and  Platt  et.al  [11]  in  that  we  studied  as  far  as  we  can  the  limit  as  e  —>  0.  In  particular  we 
computed  the  solutions  of  (2.2)  with  forcing  wavenumbers  kj  =  10,25,50,100  (corresponds 
to  £  =  0.1,0.04,0.02,0.01).  In  these  cases  the  solutions  exhibit  a  clear  separation  of  scales 
for  some  time. 

Here  we  report  two  examples  of  our  computational  results.  We  have  experimented  with 
other  sets  of  initial  data.  In  all  cases  we  observe  the  same  kind  of  phenomena  as  we  report 
below. 

In  the  first  example  we  take  the  initial  vorticity  to  be 


(2.3) 


9uq  dVQ  ,  o  X]  ■  2  2 

Wo  =  ^-^7  =  sin  Tsm  T-* 


We  computed  the  energy  spectrum  using  the  following  formula: 


(2.4)  E(k)=  £  (\u(kuk2)\2  +  \i(kuk2)\2) 

lvAi+fcJ-fcl<  I 

where  w(&i,&2)  is  the  (A'i,&2)-th  Fourier  mode  of  u.  Figure  1  is  the  log-log  plot  the  energy 
spectrum  at  times  tj  =  0.1, 0.2,  •  ■  -,2.5.  Here  £  =  0.04,  and  the  computation  is  done  on  a 
4002  grid.  It  is  clear  that  the  energy  spectrum  consists  of  two  parts:  The  high  frequency 
part  is  dominated  by  the  forcing  frequency  and  its  higher  harmonics.  The  low  frequency 
part  corresponds  to  the  large  scale  structure  of  the  flow.  Although  it  is  not  shown  in  Figure 
1,  the  energy  spectrum  at  t  -  0  decays  much  faster  than  those  plotted  there.  At  early  times 
the  low  frequency  spectrum  is  very  well  separated  from  the  high  frequency  part.  At  about 
t  =  1,  the  high  end  of  the  low  frequency  spectrum  begins  to  merge  with  the  low  end  of  the 
high  frequency  spectrum.  At  the  same  time,  the  low  end  of  t'.e  low  frequency  spectrum 


forms  an  envelope.  This  envelope  is  compared  to  a  straight  line  with  slope  equal  to  -4. 
They  are  close. 

In  the  second  example  we  choose  a  nonsmooth  initial  vorticity: 

(2.5)  w(ii,x2)  =  H(Xi)H(xi)  -  7r4, 

where  H(x)  is  the  hat  function  on  [0,27r]  with  height  n.  The  other  parameters  are  the  same 
as  in  the  first  example.  The  energy  spectrum  for  t  =  0.1, 0.2, •  •  •,  1  are  plotted  in  Figure  2. 
We  observe  basically  the  same  phenomena  as  in  Figure  1  except  that  here  the  initial  data  is 
much  less  smooth,  therefore  the  low  frequency  spectrum  settles  down  to  an  envelope  much 
faster.  The  envelope  is  also  compared  to  a  straight  line  of  slope  -4  and  is  seen  to  be  very 
close  to  that. 

Although  there  exists  a  large  literature  on  two-dimensional  turbulence  theory  (for  a 
review,  see  [2]),  the  discussions  are  much  less  sound  than  Kolmogorov’s  theory  for  three- 
dimensional  turbulence.  We  therefore  seek  our  own  explanation  of  this  seemingly  universal 
behavior  of  relaxation  to  a  k~ 4  law  for  the  low  frequency  spectrum.  In  the  next  section,  we 
will  derive  the  effective  equations  governing  the  evolution  of  the  large  scale  structures  and 
therefore  the  low  frequency  part  of  the  spectrum.  In  Section  4  we  will  study  the  effective 
equations  numerically  and  present  evidence  that  the  energy  spectrum  for  the  solutions  of 
the  effective  equations  evolve  to  a  k~4  decay  law,  even  with  a  smooth  initial  data.  The  time 
that  this  happens  is  close  to  the  time  when  the  envelope  forms  at  the  low  frequency  part  of 
the  spectrum  of  the  Kolmogorov  flows.  Beyond  this  time  the  original  Kolmogorov  flows  no 
longer  exhibit  separation  of  scales  and  the  effective  equations  cease  to  be  valid. 

§  3.  Effective  Equations  for  the  Large  Scales 

In  this  section,  we  derive  effective  equations  governing  the  evolution  of  large  scale  flow 
quantities.  Enormous  efforts  have  been  devoted  to  the  the  derivation  of  similar  equations, 
the  so-called  “Reynolds  equation”,  for  general  turbulent  flows.  The  fact  that  turbulent  flow 
has  a  continuum  range  of  scales  presents  a  serious  difficulty  to  this  objective.  Our  problem 
is  drastically  simplified  since  there  is  a  separation  of  the  scales  which  are  active  in  the  flow. 
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As  we  pointed  out  in  the  last  section,  eventually  separation  of  scales  no  longer  holds  and 
the  flow  exhibit0  a  transition  to  a  continuum  of  scales.  This  signals  the  breakdown  of  the 
effective  equation,  an  issue  which  will  be  looked  at  more  carefully  in  the  next  section. 

We  begin  with  the  following  proposition: 


u‘(x , t ) 
ve(x,  t) 


u(x,  t) 
v(x,  t) 


+  |  w(x,t ,f*)  \  +£( 


vl(xt  ~g^) 


p€{x,t)  =  p{x,t)  +  n  (x,t,  ~)  +  epi(x,t,  y)  +  ••• 

where  w(x,t,y2),  U(x,t,r/2),  ui(x,tfy2),  Pi(x,t,y2 )  etc.  are  periodic  functions  of  y2,  and 
<  w(x,  <,•)>=  0,  <  z(x,  t,  •)  >=  0,  <  II(:r,  t,  ■)  >=  0  for  every  x,  t  G  R2  x  R+ .  Substituting 
(3.1)  into  (2.2),  and  collecting  equal  orders  of  £,  we  get  a  hierarchy  of  equations.  The  0(e~l) 
equation  are: 


(3.2)  Vy 


(u  +  in)2  ( u  +  w)(v  +  z) 


+  vvn  —  Ar 


sin  3/2 


\  (u  +  w)(u  +  z)  (v+w)2  )  \  z  ) 

Vy  ■  ( W,Z)T  =  0. 

The  solutions  of  (3.2)  can  be  obtained  explicitly: 

(3.3)  w  =  -= - [sin  y2  ~  t’cos  y2]i  z  =  0,  11  =  0. 

v*  T  1 

We  next  come  to  the  0(1)  equations: 


)♦ 

+  Vt-  [ 

/ 

t 

*  / 

,  \ 

u(x, t) 
v(x,t) 


2(u  +  w)u\ 

(u  +  w)v  1  +  (v  +  z)u\ 


(u  +  w)2  (u  +  w)(v  +  z) 
(u  +  w)(v  +  z)  (v  +  w)2 


(u  +  w)v\  +  (v  +  z)ui 
2(v  +  z)v  I 


+  vxp  +  vxii  +  v  yp\ 


W  /  U 1 

=  2(VX-Vy)  +AJ 

2  /  \  *'i 


Vi  •  (u  +  w,i>  +  z)T  +  Vy  ■  (U\,l'i)T  =  0. 
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Averaging  both  sides  of  (3.4)  with  respect  to  y,  we  get  (written  in  conservative  form): 


u(x,t )  u2  uv 

+  Vi  •  + 

(3.5)  <  y  v(x,t)  y  y  uv  v 2  ) 

Vx  •  it  =  0 

Substituting  in  (3.5)  the  result  of  (3.3),  we  get 

iit  +  +  vuX2  +  (^ 

(3.6) 

4  +  Mir,  +  VVt +  ir. 


+  Vr  • 


<  w  >  <  wz  > 


<  wz  >  <  z*  > 


+  Vxp  =  0 


Ut  +  UUXl  +  VUX2  +  ( l"))x»  +  Pari  ~  0 


Vt  +  UVXi  +  vvX2  +pX2  =0 
til,  +  =  0 

This  is  the  effective  equation  governing  the  evolution  of  the  large  scale  quantities  for  the 
Kolmogorov  flow.  For  the  self-consistency  of  this  asymptotic  argument,  we  have  to  produce 
a  candidate  (ui,pi)  which  solves  the  remainder  of  (3.4).  To  this  end  we  write  (3.4)  as 

2(u  +  u;)ui  (u  +  w)vi  +  (v  +  z)u\  \  (  ux 

(3.7)  V,  •  +  VyPl  -  A  J  =  / 

^  (u  +  w)vi  +  (v  +  z)ux  2(v  +  z)vi  J  y  Vi  ) 

(3.8)  Vy  ■  (ui,  tq)T  =  -  V*  •  (u  +  w,  v  +  z)T  =  g 


(3.7)  V,. 

Ac 

(3.8) 

where 

/ 

-/  = 

\ 

-fi 

—  I 2 

( u  +  w)2  ( u  +  w)(v  +  z ) 

(u  +  in)(u  +  z)  (i  +  w)2 


u(x,  t) 
v(x,t) 


+vxp  +  vxii-2(vI-vy)  ^  J. 

In  (3.7)  and  (3.8),  x  and  t  are  viewed  as  parameters.  Since  /  and  g  only  depend  on  j/2> 
we  seek  a  solution  of  (3.7)  and  (3.8)  which  depends  only  on  j/2-  The  constructions  above 
for  w,z,U,u,v  and  p  guarantee  that  the  averages  of  /  and  g  in  t/2  over  the  period  [0,27t] 
are  zero.  Therefore  we  solve  iq  from  (3.8).  The  first  equation  of  (3.7)  will  then  give  the 
solution  uj  whereas  the  second  equation  will  give  the  solution  pj.  It  is  easy  to  see  that  if 
we  define 


u€{x,t) 

v€{x,t) 


u(x,t) 

v(x,t) 


w(x,t,  ^) 
z{x,t,^) 


U\ (x,t,^) 

iq(x,t,  ^) 


p'(x,t)  =  p(x,t)  +  ll(x,<,— )  +  £  px(x,t,—), 

then  the  equations  (2.2)  are  satisfied  by  ue,vs,pe  with  an  error  of  order  e. 

Remark  1:  The  difference  between  equation  (3.6)  and  the  usual  2-D  Euler’s  equation 
is  the  appearance  of  the  off-diagonal  term  in  (3.6).  This  means  that  the  linearized  equation 
is  weakly  well-posed:  high  frequency  perturbations  will  grow  with  a  growth  rate  which  is 
a  polynomial  in  the  wavenumber.  This  fact  motivates  us  to  look  for  finite  time  singularity 
formation  in  the  solutions  of  (3.6)  with  smooth  initial  data. 

Remark  2:  If  we  change  the  forcing  term  in  (1.1)  to  a  different  function,  the  form 
of  the  effective  equation  will  not  change.  The  off-diagonal  term  will  usually  change  to  a 
different  function.  In  particular  if  the  viscosity  is  scaled  to  vanish  faster  than  O(s),  then 
the  added  flux  in  the  effective  equation  is  singular  in  v. 

Remark  3:  We  notice  that  (3.6)  conserves  energy: 

\jtJ  2  +  ^ dx'dx 2  =  “  /  “( 2(d+l)- ^dxidx2 
=  I  t,-'WTT)dx'dx^ =  ■/*"2(P TT)dx'dXl  =  °- 

§  4.  Singularity  Formation  in  the  Solutions  of  the  Effective  Equations 

In  this  section,  we  present  numerical  evidence  which  suggests  that  the  solutions  of  (3.6) 
with  analytic  initial  data  develop  singularities  at  finite  time.  Furthermore,  at  the  time  when 
singularity  forms,  the  energy  spectrum  associated  with  the  solutions  decays  like  k~A. 

The  numerical  results  presented  below  were  computed  using  a  Fourier-collocation  method. 
We  checked  these  results  using  both  a  second  order  centered  differencing  scheme  and  a  fourth 
order  ENO  type  scheme  [14].  We  remark  that  searching  singularities  numerically  is  a  very 
subtle  undertaking.  It  is  important  to  check  the  numerical  results  using  different  meth¬ 
ods  since  each  method  is  likely  to  have  its  own  numerical  artifact  when  it  comes  to  the 
computation  of  singularities. 

Strictly  speaking,  demonstrating  finite  time  singularity  formation  is  an  analytical  prob¬ 
lem.  The  ultimate  answer  has  to  come  analytically,  either  by  giving  specific  examples  or 
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by  establishing  a  lower  bound  for  the  life  span  of  the  smooth  solutions.  Computers  are 
limited  to  finite  arithmetic  and  finite  capacities.  A  numerical  approach  can  at  best  give 
partial  evidence,  not  the  complete  answer.  Needless  to  say  that  we  have  tried  the  analytical 
approach  and  the  problem  seems  to  be  very  hard.  For  the  moment,  presenting  numerical 
evidence  is  the  best  we  can  do  to  shed  some  light  on  our  problem. 

§  4.1.  Description  of  the  Numerical  Method 

For  the  spatial  discretization,  we  used  the  Fourier-collocation  method  [1],  Roughly 
speaking,  the  differentiation  operator  is  approximated  in  the  Fourier  space,  whereas  the 
nonlinear  operations  such  as  multiplication  are  done  in  the  physical  space.  We  use  the 
intrinsic  Cray  FFT  routines  which  considerably  enhanced  the  performance  of  the  code. 

We  observed  that  even  for  the  computation  of  smooth  solutions  of  the  standard  two 
dimensional  incompressible  Euler’s  equation,  a  certain  amount  of  filtering  is  needed  for  the 
stability  of  the  numerical  method,  and  it  is  crucial  to  add  the  filters  correctly.  A  robust  way 
of  adding  the  filters  is  to  replace  the  Fourier  multiplier  ikj  for  the  differentiation  operator 
gfj  by  ikjVdkjl),  where 

(4.1)  tp(k)  =  e  Q(7Vr)  ,  for  \k\  <  N 

Here  N  is  the  numerical  cutoff  for  the  Fourier  modes,  mj  is  the  order  of  the  filter,  a  is 
chosen  so  that  <p(N)  =  e~a  =  machine  accuracy.  The  machine  accuracy  on  Cray-YMP  with 
single  precision  is  roughly  10~14.  The  accuracy  of  such  an  approximation  scheme  depends 
on  the  parameter  mj.  For  smooth  functions  f(x),  we  have 

(4.2)  \\f\x)-DNf(x)\\  =  0(N-Tnf) 

With  DnJ  =  F~l(ikp>(\k\))F f  is  the  numerical  approximation  of  f'(x).  F  denotes  the 
Fourier  transform  operator.  Unless  otherwise  stated,  the  results  presented  in  §  4.2  were 
computed  with  mj  =  14. 

For  the  temporal  discretization,  we  used  the  third  order  Runge-Kutta  method  designed 
in  [14].  We  used  the  third  order  version  since  it  only  requires  three  auxiliary  arrays,  whereas 
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the  fourth  order  version  requires  five  auxiliary  arrays.  We  take  initial  data  that  is  periodic 
with  period  D  where  D  =  [0, 2;r]  x  [0, 2?r]. 

§  4.2.  Numerical  Results 

We  present  our  numerical  results  for  the  case  when  the  initial  vorticity  is  the  same  as 
in  the  first  example  of  Section  2:  uq  =  sin2  ^  sin2  ^  —  it2. 

As  a  first  hint  for  the  singularity  formation,  we  plotted  the  time  history  of  the  quantity 
Hi(t)  =  fD\Vijj\2dxdy.  Figure  3  suggests  that  as  IV  — ♦  +oo,  H\(t)  grows  without  bound, 
for  t  >  0.9. 

This  information  alone  cannot  be  used  as  a  solid  evidence  for  singularity  formation, 
since  it  only  says  that  the  quantity  H\(t)  is  poorly  resolved  for  t  >  0.9.  Next  we  plotted 
in  Figure  4  the  time  evolution  of  the  energy  spectrum  (on  a  log-log  scale)  associated  with 
the  numerical  solutions  obtained  with  three  different  resolutions:  1282,  2562  and  4002.  We 
observe  that  at  t  =  1,  a  bump  develops  near  the  cutoff  frequency  on  the  energy  spectrum. 
This  indicates  that  at  this  time,  the  numerical  solutions  are  not  very  well  resolved.  We  also 
observe  that  the  energy  spectrum  at  t  =  0.9  is  well  represented  by  the  numerical  solutions 
obtained  on  the  1282  and  2562  grids,  and  it  fits  very  well  with  a  straight  line  with  slope 
-4.  This  is  a  stronger  indication  that  at  this  time,  the  energy  spectrum  of  the  solution  no 
longer  decays  exponentially,  but  decays  algebraically  as  k~ 4. 

To  determine  more  precisely  when  algebraic  decay  sets  in,  we  use  the  techniques  pro¬ 
posed  by  Sulem,  Sulem  and  Frisch  [17].  Let  f{z)  be  a  function  which  is  analytic  in  the  strip 
|Inu|  <  a,  and  has  a  branch-cut  singularity  z^  on  |Imz|  =  a,  then  its  Fourier  transform 
f(k)  decays  like 

(4.3)  'f{k)~Ck-°e-ak,  |fc|-+oo. 

Imagine  the  following  scenario  for  the  singularity  formation  in  the  solutions  of  (3.6). 
Instead  of  solving  (3.6)  in  the  real  space  R 2,  we  solve  (3.6)  in  the  complex  space  (zj,  Z2)  £  C2 . 
Initially  the  solution  is  an  entire  function.  At  later  time  the  solution  develops  branch  point 
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of  order  (3(t)  at  z(t)  with  |Ini2(f)|  =  a(t)  >  0,  but  it  is  analytic  in  the  strip  1 1 m 2r |  <  a{t). 
If  this  branch  point  travels  to  the  real  axis  at  a  finite  time  t* ,  then  the  solution  of  (3.6), 
solved  in  the  real  space,  develops  a  singularity  at  t  =  t*. 

This  scenario  has  been  proposed  for  a  number  of  problem  in  fluid  mechanics,  notably 
the  formation  of  singularity  in  the  evolution  of  a  vortex  sheet,  (see  [8,  3,  13]).  In  this  case, 
there  are  strong  numerical  as  well  as  analytical  evidence  which  confirms  this  picture. 

To  check  whether  this  provides  a  plausible  picture  for  our  problem,  we  checked  numer¬ 
ically  the  validity  of  proposition  4.3  for  the  energy  spectrum.  To  do  that  we  choose  three 
consecutive  wave  numbers  (k  -  1  ,k,k  +  1),  and  solve 

(4.4)  E(k,  t)  =  Ck-Pe-a~k,  for  k  =  k  -  l,k,  k  +  1, 

with  double  precision  arithmetic  (on  Cray-YMP).  At  each  fixed  time,  we  get  three  functions 
C\k,  t),  p(k,  f),  a{k,  t).  These  functions  depend  on  k.  In  order  that  (4.3)  be  a  good  estimate 
for  the  energy  spectrum  of  the  solutions  of  (3.6),  asymptotically  as  k  — ►  oo,  the  functions 
C(k,t),  i3(k,t)  and  a(k,t)  should  not  depend  on  k  in  this  limit.  In  Figure  6  we  plotted  the 
functions  a{k,t)  from  t  =  0.4  to  t  =  1  with  constant  time  increments.  The  computation  is 
done  on  a  2562  grid.  It  is  clear  from  Figure  5  that  for  fixed  t ,  a(k,  t)  is  almost  a  constant  for 
20  <  k  <  60.  For  k  <  20,  a (k,t)  exhibits  fluctuations  since  k  is  not  yet  in  the  asymptotic 
regime.  For  80  <  k  <  128,  the  energy  density  is  at  the  order  of  10~32  which  is  below  machine 
accuracy  even  if  double  precision  is  used.  These  wavenumbers  are  needed  only  to  enforce 
smooth  decay  for  the  Fourier  coefficients  down  to  the  machine  accuracy.  Figure  6  displays 
similar  results  for  a  computation  using  a  4002  grid.  The  fluctuations  at  high  wavenumber 
is  a  manifestation  of  trying  to  fit  machine  zeros  with  proposition  4.3.  The  functions  C(k,t ) 
and  (i{k,t)  exhibit  the  same  feature. 

Now  we  can  average  the  values  of  a  over  the  intermediate  wave  numbers  and  safely 
regard  the  averaged  value  (as  a  function  of  t)  as  the  width  of  the  analyticity  strip  of  the 
solution  at  that  time.  This  function  will  be  denoted  by  a(t).  The  corresponding  function 
for  (3  will  be  denoted  by  f3(t).  Displayed  in  Figure  7  is  the  function  o(t)  for  two  different 
numerical  resolutions:  2562  and  4002.  Results  obtained  on  boMi  grids  indicate  that  the 
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width  of  the  analyticity  strip  for  the  solution  vanishes  at  t*  0.9,  suggesting  a  singularity 
formation  at  this  time.  The  corresponding  results  for  (3  are  displayed  in  Figure  8.  As  t  —*  t* , 
j3(t)  converges  to  4. 

The  other  numerical  parameter  in  our  method  is  the  order  of  the  filter.  Figures  9  and  10 
display  the  results  for  <5,  (3  as  we  vary  the  order  of  the  filter.  We  see  that  the  approximated 
values  of  a,  (3  are  not  sensitive  to  the  change  of  the  filters.  However  when  the  filter  is  too 
weak  (to/  =  15),  the  numerical  solution  deteriorates  drastically  at  a  time  before  <*,  since 
it  can  no  longer  handle  the  large  gradients  developed  in  the  vorticity  field.  The  computed 
value  fluctuates  much  more,  therefore  the  result  is  not  shown  in  Figure  10.  This  is  not 
unexpected  since  in  general  (3  is  much  more  sensitive  than  a. 

To  get  an  idea  about  the  nature  of  the  singularity,  we  display  in  Figures  11  and  12  the 
time  evolution  of  u  and  u  at  y  =  4^-  with  t  =  0.3, 0.6, 0.9.  At  t  =  0.9,  a  cusp  seems  to  have 
formed  on  the  profile  of  ui.  Contour  plots  of  u  show  that  these  cuspy  structures  occur  on 
isolated  lines. 

The  behavior  of  the  solutions  reported  above  seems  to  be  generic.  It  occurs  in  other 
calculations  we  did  using  different  initial  data. 

§  5.  Concluding  Remarks 

As  we  mentioned  in  the  introduction,  the  Kolmogorov  flow  belongs  to  a  large  class  of 
problems  which  exhibit  large  scale  instability.  A  typical  example  in  three  dimension  is 
the  ABC  flow  which  is  known  to  be  unstable  to  large  scale  perturbations  [7].  In  principle 
the  ideas  and  methods  presented  above  should  apply  to  these  problems  also,  although  to 
actually  carry  out  the  program  is  difficult  mainly  because  the  computational  cost  increases 
drastically  for  three-dimensional  problems.  However,  it  is  clear  that  such  instabilities  will 
trigger  the  transport  of  energy  from  small  scales  to  large  scales,  thereby  effecting  the  mean 
flow  quantities. 

One  can  also  attempt  to  study  the  interaction  of  different  scales  by  studying  the  solutions 
with  two-scale  initial  data.  Formally  one  can  use  multiple  scale  asymptotics  to  derive 
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effective  equations  for  the  large  scales.  This  is  done  in  [5,  4].  However,  a  closer  look  at 
these  effective  equations  reveals  that  the  linearized  equations  are  exponentially  ill-posed: 
high  frequency  perturbations  grow  at  a  rate  that  is  exponential  in  the  wavenumber.  The 
reason  for  this  is  very  simple.  The  incompressibility  condition  forces  the  microstructure 
to  be  shear  flows.  In  the  limit  of  infinite  scale  separation,  the  microstructure  appears  like 
stacks  of  vortex  sheets  when  viewed  at  macro-scales.  The  linear  ill-posedness  is  simply  a 
manifestation  of  the  Kelvin-Helmholtz  instability. 

In  such  situations  multiple  scale  asymptotics  does  not  provide  useful  information  since 
it  is  based  on  the  hypothesis  of  scale  separation.  For  the  Kolmogorov  flow,  we  are  saved 
since  there  is  a  fixed  time  interval  in  which  the  scales  are  separated  and  we  can  get  useful 
information  from  two-scale  expansions.  As  one  might  expect,  beyond  the  critical  time  t*,  the 
two  scales  in  the  Kolmogorov  flow  merge  and  the  effective  equation  cease  to  be  valid.  This  is 
checked  numerically  with  the  following  procedure.  The  effective  equation  (3.6)  is  integrated 
numerically  using  a  finite  difference  method  with  a  built-in  nonlinear  numerical  viscosity 
to  avoid  oscillations  coming  from  large  vorticity  gradient  [14].  The  energy  spectrum  of  the 
numerical  solution  is  found  to  be  decaying  increasingly  slower.  This  differs  sharply  from  the 
energy  spectrum  of  the  original  Kolmogorov  flow  since  the  latter  approaches  an  asymptote 
after  the  critical  time. 
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Figure  1.  Log-log  plot  of  the  time  evolution  of  the  energy  spectrum  with  initial  data  (2.3) 
at  i  =  0.1, 0.2,  •  ■  •,  2.5.  £  =  0.04.  t  increases  upward.  The  energy  spectrum  is  compared  to 
a  line  of  slope  -4. 


at  t  =  0.1, 0.2,  ■■■,!.£  =  0.04.  t  increases  upward.  The  energy  spectrum  is  compared  to  a 
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Figure  3.  From  bottom  to  top,  time  history  of  Hi(t )  on  computational  grids:1282, 2562, 400 


Figure  4.  Log-log  plot  of  the  time  evolution  of  the  energy  spectrum  computed  on  grids: 
1282,2562and  4002.  From  bottom  to  top:  t  =  0.1, 0.2, 1. The  energy  spectrum  is  com¬ 
pared  to  a  line  of  slope  -4. 
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Figure  10.  f3(t)  computed  on  a  2562  grid  with  different  orders  of  filter.  The  line  with  circles 
is  computed  with  mj  =  5.  The  other  line  is  computed  with  mj  =  10.  The  result  for 
m/  =  15  fluctuates  and  is  not  shown  here. 
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